####Please ensure that all of the below libraries are installed on replicating machine####

library(maptools)
library(spdep)
library(rgdal)
library(maps)
library(foreign)
library(rgeos)
library(gdata)
library(McSpatial)
library(biganalytics)
library(bigmemory)
library(Matrix)

####Replication for Figure 2 in the main paper####

conf=.95
title=""
xlabel="Region Contains Insurgency"
ylabel="Incidence Rate Ratio of Attacks"
factor_labels=c("No","Yes")


# Get coefficients of variables
beta_1 = -.2521635
beta_3 = .0545515

# Create list of moderator values at which marginal effect is evaluated
x_2 <- c(0,1)

# Compute marginal effects
delta_1 = exp(beta_1 + beta_3*x_2)

# Compute variances
var_1 =   .00784781  + (x_2^2)*-.00043512   + 2*x_2*.01686156  

# Standard errors
se_1 = se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window
plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")

# Plot points of estimated effects
points(x=x_2, y=delta_1, pch=16)

# Plot lines of confidence intervals
lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")

# Label the axis
axis(side=1, at=c(0,1), labels=factor_labels)

# Add a dashed horizontal line for zero
abline(h=1, lty=3)


####Replication for Figure 1 in the appendix####

###Note: For proper replication, ~PATH for each of the shapefiles should be changed
###to the actual location of the shapefiles on the user's computer

###Adjacency Matrix:Algeria###
setwd("~PATH/DZA_adm")
map<-readOGR(dsn = ".", layer='DZA_adm1') 
plot(map)
map$NAME_1
dza.nb = poly2nb(map)
adj1 = nb2mat(dza.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
dza.adj<-adj
dza.inv<-ginv(dza.adj)


###Adjacency Matrix:Nigeria###
setwd("~PATH/NGA_adm")
map<-readOGR(dsn = ".", layer='NGA_adm1') 
plot(map)
nga.nb = poly2nb(map)
adj1 = nb2mat(nga.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
nga.adj<-adj
nga.inv<-ginv(nga.adj)




###Adjacency Matrix:Nepal###
setwd("~PATH/NPL_adm")
map<-readOGR(dsn = ".", layer='NPL_adm2') 
plot(map)
npl.nb = poly2nb(map)
adj1 = nb2mat(npl.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
npl.adj<-adj
npl.inv<-ginv(npl.adj)


###Adjacency Matrix:Peru###
setwd("~PATH/PER_adm")
map<-readOGR(dsn = ".", layer='PER_adm1') 
plot(map)
per.nb = poly2nb(map)
adj1 = nb2mat(per.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
per.adj<-adj
per.inv<-ginv(per.adj)

###Adjacency Matrix:Syria###
setwd("~PATH/SYR_adm")
map<-readOGR(dsn = ".", layer='SYR_adm1') 
plot(map)
syr.nb = poly2nb(map)
adj1 = nb2mat(syr.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
syr.adj<-adj
syr.inv<-ginv(syr.adj)




###Adjacency Matrix:Turkey###
setwd("~PATH/TUR_adm")
map<-readOGR(dsn = ".", layer='TUR_adm1') 
plot(map)
tur.nb = poly2nb(map)
adj1 = nb2mat(tur.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
tur.adj<-adj
tur.inv<-ginv(tur.adj)


###Adjacency Matrix:Uganda###
setwd("~PATH/UGA_adm")
map<-readOGR(dsn = ".", layer='UGA_adm1') 
plot(map)
uga.nb = poly2nb(map)
adj1 = nb2mat(uga.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
uga.adj<-adj
uga.inv<-ginv(uga.adj)



###Adjacency Matrix:Yemen###
setwd("~PATH/YEM_adm")
map<-readOGR(dsn = ".", layer='YEM_adm1') 
plot(map)
yem.nb = poly2nb(map)
adj1 = nb2mat(yem.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
yem.adj<-adj
yem.inv<-ginv(yem.adj)




###Adjacency Matrix:Thailand###
setwd("C:/Users/konst/Desktop/Inequality/THA_adm")
map<-readOGR(dsn = ".", layer='THA_adm1') 
plot(map)
tha.nb = poly2nb(map)
adj1 = nb2mat(tha.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
tha.adj<-adj
tha.inv<-ginv(tha.adj)


###Adjacency Matrix:Russia###
setwd("~PATH/RUS_adm")
map<-readOGR(dsn = ".", layer='RUS_adm1') 
plot(map)
rus.nb = poly2nb(map)
adj1 = nb2mat(rus.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
rus.adj<-adj
rus.inv<-ginv(rus.adj)


###Adjacency Matrix:Iraq###
setwd("C:/Users/konst/Desktop/Inequality/IRQ_adm")
map<-readOGR(dsn = ".", layer='IRQ_adm1') 
plot(map)
irq.nb = poly2nb(map)
adj1 = nb2mat(irq.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
irq.adj<-adj
irq.inv<-ginv(irq.adj)


###Adjacency Matrix:Cambodia###
setwd("~PATH/KHM_adm")
map<-readOGR(dsn = ".", layer='KHM_adm1') 
plot(map)
map$NAME_1
khm.nb = poly2nb(map)
adj1 = nb2mat(khm.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
khm.adj<-adj
khm.inv<-ginv(khm.adj)


###Adjacency Matrix:Sri Lanka###
setwd("~PATH/LKA_adm")
map<-readOGR(dsn = ".", layer='LKA_adm1') 
plot(map)
lka.nb = poly2nb(map)
adj1 = nb2mat(lka.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
lka.adj<-adj
lka.inv<-ginv(lka.adj)



###Adjacency Matrix:Myanmar###
setwd("~PATH/MMR_adm")
map<-readOGR(dsn = ".", layer='MMR_adm1') 
plot(map)
mmr.nb = poly2nb(map)
adj1 = nb2mat(mmr.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
mmr.adj<-adj
mmr.inv<-ginv(mmr.adj)



###Adjacency Matrix:Pakistan###
setwd("~PATH/PAK_adm")
map<-readOGR(dsn = ".", layer='PAK_adm1') 
plot(map)
pak.nb = poly2nb(map)
adj1 = nb2mat(pak.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
pak.adj<-adj
pak.inv<-ginv(pak.adj)



###Adjacency Matrix:Afghanistan###
setwd("~PATH/AFG_adm")
map<-readOGR(dsn = ".", layer='AFG_adm1') 
plot(map)
afg.nb = poly2nb(map)
adj1 = nb2mat(afg.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
afg.adj<-adj
afg.inv<-ginv(afg.adj)



###Adjacency Matrix:Angola###
setwd("~PATH/AGO_adm")
map<-readOGR(dsn = ".", layer='AGO_adm1') 
plot(map)
ago.nb = poly2nb(map)
adj1 = nb2mat(ago.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
ago.adj<-adj
ago.inv<-ginv(ago.adj)

###Adjacency Matrix:Burundi###
setwd("~PATH/BDI_adm")
map<-readOGR(dsn = ".", layer='BDI_adm1') 
plot(map)
bdi.nb = poly2nb(map)
adj1 = nb2mat(bdi.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
bdi.adj<-adj
bdi.inv<-ginv(bdi.adj)



###Adjacency Matrix:India###
setwd("~PATH/india")
map<-readOGR(dsn = ".", layer='IND_adm1') 
plot(map)
ind.nb = poly2nb(map)
adj1 = nb2mat(ind.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
ind.adj<-adj
ind.inv<-ginv(ind.adj)




###Adjacency Matrix:Colombia###
setwd("~PATH/colombia")
map<-readOGR(dsn = ".", layer='COL_adm1') 
plot(map)
col.nb = poly2nb(map)
adj1 = nb2mat(col.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
col.adj<-adj
col.inv<-ginv(col.adj)




###Adjacency Matrix:Philippines###
setwd("~PATH/phillippines")

map<-readOGR(dsn = ".", layer='Regions') 
plot(map)
phl.nb = poly2nb(map)
adj1 = nb2mat(phl.nb, style="W",zero.policy=T)
adj = as.matrix(bdiag(adj1,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))  
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
adj = as.matrix(bdiag(adj,adj1))
phl.adj<-adj
phl.inv<-ginv(phl.adj)

####Each adjacency matrix is compiled and replicated twenty times, to create a sufficient
####amount of matrices to reflect the time span of the data. The matrices can then be put 
####together using the below code.

####Putting weight matrices together####
spw<-NULL
spw<-bdiag(afg.adj,afg.adj,afg.adj)
spw<-bdiag(spw,ago.adj,ago.adj,bdi.adj)
spw<-bdiag(spw,col.adj,col.adj,dza.adj)
spw<-bdiag(spw,ind.adj,ind.adj,ind.adj,ind.adj,ind.adj,ind.adj,ind.adj,ind.adj,ind.adj,ind.adj)
spw<-bdiag(spw,irq.adj,khm.adj,lka.adj,mmr.adj,nga.adj,nga.adj)
spw<-bdiag(spw,npl.adj,npl.adj,pak.adj,pak.adj,pak.adj,pak.adj)
spw<-bdiag(spw,per.adj,per.adj,phl.adj,phl.adj,phl.adj,rus.adj)
spw<-bdiag(spw,syr.adj,syr.adj,tha.adj,tur.adj,tur.adj,uga.adj,uga.adj)
spw<-bdiag(spw,yem.adj,yem.adj,yem.adj)
spw<-as.matrix(spw)
listw<-mat2listw(spw)
rm(listw)
spw.sum<-colSums(spw, na.rm = FALSE, dims = 1)
s.lag<-

###Moran's I Test for Spatial Dependence###
setwd("~PATH")
# Note: the below data is a balanced form of the panel data in the main file, it has been cleaned up and some observations have been deleted.
main_data<-read.dta("ii_rebelterr_balanced.dta")
library(MASS)
MoransI <-moran.test(main_data$iattacks_t, listw=listw,zero.policy = TRUE)
MoransI
moran.plot(main_data$iattacks_t, listw=listw,zero.policy = TRUE,xlab = "Terrorist Attacks in Region", ylab = "Spatially Lagged DV")
##Calculate spatially lagged coefficient###
terror<-main_data$iattacks_t
terror<-t(terror)
s.lag<- as.vector(crossprod(spw,terror))
summary(s.lag)
main_data.new<-cbind(main_data,s.lag)
write.dta(main_data.new,'ii_rebelterr_balanced2.dta')
# The resultant .dta file should produce a data-set with a spatially lagged DV that then runs using STATA code in the .do file under Table 13 of the appendix.
